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ABSTRACT 


The motion in a ternary system consisting of a central body, A, close companion B, and distant 
companion, c, is considered and the perturbations of B due to C are determined. The osculating 
orbit of B relative to A can be an eccentric ellipse, not extremely elongated (say, e £ 0.3), and the 
orbit of c around the center of mass of A and Bis a fixed ellipse which might have a large ec- 
centricity, e' . The mutual inclination, I of the orbits of B and c can also be large. The possibility 
of solving the problem without limitation to small e' was indicated earlier by Brown. His basic 
small parameter is K = [^(iiq +m 2 +m 3 )J (n'/n) (l - e ' 2 ) _3/2 . The use of A. and the introduction of 
the true anomaly f' of C as the independent variable permit one to take the effect of e' into account 
in a closed form without developing series in powers of e' . 

The solution pursued here is a development of perturbations into trigonometric series with 
arguments linear in f ' and similar in form to the standard arguments of the lunar theory. To this 
end we consider the squares of k 9 of the parallax, and of the mean eccentricity of B as small 
parameters. 

The problem is solved by applying the Krylov -Bogolubov method to eliminate the significant 
short-period effects from the Milankovich differential equations for the general perturbations of 
the Laplacian vector e and the areal velocity c. Thus, the elements which are affected only by the 
long period perturbations must then be established. 

We conclude that for small eccentricity there are resonances at I ^ 39° and at I ~ 141°, and 
that the proposed trigonometric series solution cannot be done outside the intervals 0° ^ l 0 <^39° 
and 180° = l 0 >^141°. This fact gives rise to a discussion of the convergence and limits of applica- 
bility of the Delaunay series. Some of our results are compared with those of Delaunay. We con- 
clude that the Krylov -Bogolubov method can improve the Delaunay theory because it leads to more 
compact series and their convergence can be investigated more easily. In addition, they apply to 
any orbital inclination up to the critical value. 
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APPLICATION OF THE KRYLOV-BOGOLUBOV 
METHOD TO SOLUTION OF THE STELLAR THREE-BODY PROBLEM 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 


The motion in a ternary system consisting of the central body, A, the close companion B, and 
the distant companion, C, is considered in the present paper. The masses of the components are 
m i , m 2 , and m 3 , respectively. 

We assume that the osculating orbit of B relative to A is an eccentric ellipse, not extremely 
elongated (say, e $ 0.3), and that the orbit of C around the center of masses of two close companions 
A and B is a fixed ellipse which might have a large eccentricity, e\ The mutual inclination I of the 
orbits of B and C can also be large. In the further exposition we shall discuss a limitation which 
must be imposed on I in order to avoid a resonance condition. 


We have to determine the perturbations of B as caused by C. The problem so formulated is 
called the stellar three-body problem because the conditions described above exist in some triple 
star systems. The case of the lunar satellite, because neither e nor I must be small, evidently 
represents a special case of the stellar problem. The role of C is played by either the Earth or 
the Sim. If e, e' , and I are small we have the classic lunar problem; thus the importance of the 
general stellar problem is evident, and attempts should be made to solve it in an analytic or semi- 
analytic way. 


The works of Kozai and Kovalevsky (References 1 and 2, respectively) on the problem of mo- 
tion of the lunar satellite deserve to be mentioned. They solved the problem by assuming that 
e' =0. The possibility of solving the problem without this limitation was indicated by Brown in 
his work on the stellar three-body problem (Reference 3). The basic small parameter in Brown's 
theory is 


k = 


IT (l-e‘ 2 )- 3 ' 2 , 


where n ' is the mean motion of c and n is the mean motion of B . In the theory of Delaunay (Refer- 
ence 4) the basic parameter is n'/n, but the use of k, and also the introduction of the true anomaly 
f ' of C as the indendent variable, permit one to take the effect of e' into account in a closed form 
without resorting to the development of series in powers of e' . 
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The form of the solution we are pursuing here is the development of perturbations into trig- 
onometric series with the arguments linear in f # and similar in form to the standard arguments 
of the lunar theory. In order to achieve such a solution we consider the square of the mean eccen- 
tricity of B, k, and the parallax as small parameters. This limitation is justified by the circum- 
stance that the orbits of lunar satellites with large eccentricities will be highly unstable. In addi- 
tion the fact that the Moon is a body and not a point, together with the assumption that the parallax 
is small, imposes a certain upper limit on the eccentricity of B. The mutual inclination I and e # , 
however, are not presupposed to be small. 

In keeping the orbit stable, the long-period perturbations are more important than the short- 
period ones. This fact is the foundation of every satellite theory which makes use of the method of 
variation of astronomical constants. The way to solve the problem is to eliminate gradually all of 
the significant short-period perturbations from the differential equations. Thus the elements which 
are affected only by the long- period perturbations must be established. The establishment of such 
elements and the formation of the differential equations governing their variations require the solu- 
tion of a chain of linear partial differential equations of the first order. After the differential equa- 
tions for the new elements are finally formed, we can integrate them using the Poincare' small pa- 
rameter method (Reference 5), provided the eccentricity is not large. 

The elimination of the short-period terms is facilitated if in the original differential equations 
the coefficients of the derivatives of the disturbing function are linear in the elements. For this 
reason the Laplacian vector e and the areal velocity c represent a convenient choice. Their gen- 
eral perturbations are governed by the differential equations of Milankovich (Reference 6). The 
vectorial products contained in these equations are the linear functions of e and c. The inconven- 
ience caused by the division by e in the equation for de/dt disappears after the elimination of the 
short-period effects. These elements have a direct kinematic meaning and with their help we can 
easily visualize the motion of the orbit in space. They are not canonical and consequently it is 
convenient to apply the method of Krylov-Bogolubov (Reference 7) to solve the problem. 

The differential equations for the elements affected only by the long-period perturbations can 
be written as a quasi-linear system. We can solve it either by developing the solutions into power 
series in a small parameter or by applying the method of successive approximations. In recent 
years we have applied the digital computer to the development of the semianalytic lunar theories 
of the Hansen type (Reference 8). The machine develops the disturbing function, forms the differ- 
ential equations, and integrates them by the process of successive approximations. In this way 
we produced the theories of lunar satellites and of the xth satellite of Jupiter. As a consequence 
of this experiment it is suggested that one substitutes the numerical values of the constants of in- 
tegration at the outset and then let the machine obtain the solution in the form of trigonometric 
series with purely numerical coefficients by means of successive approximations. At each ap- 
proximation we shall obtain c x , the mean motion of the longitude of the perigee, from a quadratic 
equation. 
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If the mean eccentricity e is zero, then, as we shall see 

c x = | * -i + 0(K 2 e 2 ) , 

where k is the sine of the mean inclination i 0 and e is the mean eccentricity. From this we conclude 
that for small eccentricity there is a resonance at l 0 ~ 39° and also at l 0 ~ 141°. We also conclude 
that the development of perturbations of e and c into trigonometric series with the standard argu- 
ments of the lunar theory cannot be done outside the intervals 0°= I 0 < ~39° and 180° 1 1 0 > ~141°. The 
existence of this critical inclination in the special case e' =0 was pointed out by Kevorkian (Ref- 
erence 9) and Kozai (Reference 10). The fact that it exists also for e' / 0 gives rise to the discus- 
sion about the convergence and limits of applicability of the series of Delaunay (Reference 4). In a 
further exposition we shall compare some of our results with the results of Delaunay. 


BASIC DIFFERENTIAL EQUATIONS 


Let P be the unit vector directed toward the osculating pericenter of B and R be the unit vector 
normal to the osculating orbit plane of B. For the components of P and R relative to an inertial 
system located in the orbit plane of C, we have a set of standard formulas: 


P 1 = cos 2 


1 

2 


cos 77 + sin 2 


1 

2 


cos (77 - 211), 




sin 77 - sin^ 


1 

2 sin (77 - 20), 


> 


(i) 


P sin (n - Q) sin I ; 

J 


Rj - sin O sin I , 


R 2 = - cos 0 sin I , 


y ( 2 ) 


R 3 = cos I. 


Then 


e - eP * 


J 


(3) 


c = 1 - e 2 ) R 


(4) 
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If the equation of motion of the satellite is taken in the normalized form- 


d 2 r 
dt 2 

then the equations of Milankovich have the form 


— + grad R, 


(5) 


de — BR 1 - e 2 BR c 2 — BR 
-ex + c x + n —— e 


dt 


Be 


Be 


dt 


( 6 ) 


dc ^ BR _ BR 

= c x + e x _ ^ 

d t Be Be 


and 


dt 


BR 


(7) 


(8) 


Decomposing e, c, BR/ 3 e,and 3 R/ 3 c along the axes of the inertial system, with the x- and y-axes 
in the orbit plane of C, we have the system of Milankovich' s scalar equations: 


de 


A = e 


dt 


BR 


BR 


2 ^ " 3 rr | + 

Bc 3 3c 2 


Lzll (c — -c 3 ^L\ + n 

c 2 


BR BR \ . _ c 2 BR 

e 2 B-t *’ 


(9) 


de„ 


- e 2 

dt 


/ BR 

BR \ 

1 - e 2 
+ 

c 2 

/ 3 R 

BR \ 

(e 3R 

BR \ 

1 1 - (r ** 

BR 

_ c 


2 3c iJ 

1 2 
c z 

V 1 3e 2 

2 ^ 


B-t 


(10) 


( 11 ) 


dc 

dt 

dc 2 

dt 


I = c 


BR BR \ / BR 

3c 3 3 ac 2 / \ 3e 3 


BR 


BR BR \ / BR BR 


( 12 ) 


(13) 


£■(• 


BR 

dc n 


BR 

Be! 


+ e. 


BR 

bT 


BR 

Be, 


' 1 / \ ^2 

We shall make use of the complex elements 

«! = e x + ie 2 > h i = " C 2 + ic i’ 


g 2 = e x - ie a , h 2 = -c 2 - ic^ 


Ba = - 1 e v 


h 3 C 3 * 


(i=/T) 


(14) 
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instead of the real elements e 1? e 2 , e 3? c 1? c 2 , and c 3 . Their introduction permits us to operate 
with power series instead of trigonometric series and to put the result in compact and symmetri- 
cal form. There is an affinity between the use of the complex coordinates in HilFs lunar theory 
and the use of the complex vectors in the present theory. From Equations 1 and 2 we have 


g 1 = e cos 2 -L exp i 77 + e sin 2 L exp i (- n + 2 Q) 


- e(P l + i P 2 ) , 


(15) 


cos 2 — exp(- i77) + e si n 2 — exp i (tt - 2 12 ) 


= e(P 1 - iP 2 ) 


g 3 = — e sin I [ e xp i (- 77 + 12 ) - e xp - i (- 77 + 12)] 


i eP 3 , 


hj - \/a(l - e 2 ) sin I exp i fl - \/a(l - e 2 ) (- R 2 + iRj), 


and 


h 2 = \fa(l - e 2 ) sin I exp(-(l) = \fa(l - e 2 ) (- R 2 - iRj), 


h 3 = ya(l - e 2 ) c o s I = \/a( 1 - e 2 ) R 3 . 

We have from Equations 15 to 20: 


3£ _ 3R _.BR 
3gj 3e t 3e 2 


_ 3R_ = 3R +i 3R_ 
3e. 


3g 2 3e t 


3R = + . 3R 

a S 3 ^ 


3R = _ _3R^ _ . 3R 
3hj 3 c 2 3c 3 


2 3R = 3R . _3R_ 
^2 = " ^ + 1 


(16) 

(17) 

(18) 

(19) 

( 20 ) 

( 21 ) 

(22) 

(23) 

(24) 

(25) 
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and 


(26) 


BR = BR 

3h 3 Be, 

Equations 9 to 14 in the new variables are: 


s L — - g, ^ (2 h 3 I 1 4. hl ~) + n 2* gl , (27) 

dr \ 83 3h, 81 ah 3 / c 2 V 3 3g 2 1 ^ 3 / e 2 BA. 


dgj 

dr 


BR BR \ 1 - e 


= 2g 3 ^ + g. 


Bhj 3h 3 


, , BR . BR \ n c 2 3R o 

3 3g t 2 3g J e 2 3\ 2 


(28) 


dgs 

dr 


BR _ BR \ 1 - e 2 ( u BR tU BR \ „ c 2 BR 

3 «i 


g, + g, — 

1 Bhj Bh 


1 — -"2 371 T " * 3X ° 3 ’ 


(29) 


dh 


1 = 2h 


and 


where we put 


dr 


dr 


dh. 


BR 

^h 2 

- hl 

' dh 3j 



3R BR \ 

+ g, ’ 

3g 2 3g J 

(30) 

BR 

B hl 

2 

-I 

( 

2g 3 

< 

!L +§2 2L\. 

3g 2 3gJ 

(31) 


= (h, ^ - h, 


dr l 1 Bh 


BR 

Bh, 


+ — - S; 


BR 


1 2 ** 2 / 


T = it. 

X = il. 

This system of differential equations admits a first integral 

gj h, - g 2 h t + 2g 3 h 3 = 0, 

which is equivalent to the statement that 

P • R = 0. 


(32) 

(33) 

(34) 

(35) 


If we eliminate from R the short- period terms containing l in the argument, then the new 
disturbing function [R] will contain only V, n, and fi in the argument. Thus [R] will yield only 
the long- period and secular effects in the elements. 

The elimination of the short-period terms can be performed either by the method of Von Zeipel 
(Reference 11) or by the Krylov- Bogolubov method. In the treatment of the perturbations of the 
noncanonical elements the Krylov- Bogolubov method should be used. 

In the problem of motion of a close satellite we can set, with an accuracy up to X 3 , 


[ r] = — f Kd-e, 
2n Jo 


( 36 ) 
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where the long-period terms having A 4 as a factor are neglected. With the elimination of the ef- 
fects depending upon Z the element a becomes invariable. Then it is more convenient to define 
h x , h 2 , and h 3 by 

hj = v/l - e 2 (- R 2 + i Rj), 


h 2 = v/l - e 2 ( - Rj - i Rj)> 


h 3 =/TTe2 r 


(37) 

(38) 

(39) 


than by Equations 18 to 20. 


Now we set 


c = v/l - e 2 R. 

With these new definitions Equations 27 to 32 take the form: 


(4') 


dgj 

dr 


30 


30 


l h i 


an 


an 




(40) 


dg 2 _ o / 3n h 3 H\ /, an an 

dr “ ( Bs ah t 3 B g J ( 2 Bg 3 + g; 


3h, 


(41) 


dg 3 / Bn Bn 

= 6i + g 2 ’ 

dr V 1 3h 2 


.an . Bn 

h, + h 

1 Bg t 3g 2 


(42) 


dh 


1 = 2 Ik — - g. (h 


dr \ ' 3 Bh, 3 3g. 


an an 

i — +8, — 

2/ \ 3h 3 ^ g 3 


(43) 


dh 2 „ /. Bn \ (u 

— = - 2 ( h, _ + g, — ) + I h, — - g 2 


dr 


Bhj 3 8i 


ah, 3g, 


(44) 


and 


where we put 


^3 = / an 

dr \ 1 3g t 


an \ . an . an 

62 W + V 1 ^ " 2 ^ 


(45) 


From 


and 


0 = 


[R] 

/S 


e 2 = gj g 2 - g 3 2 


1 - e 2 = hj h 2 + h 3 2 , 


(46) 

(47) 

(48) 
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r 


we conclude that in addition to the integral 33 the system of Equations 40 to 45 admits also the 
integral 

Bi g 2 + h l h 2 - 83 + h 3 2 = 1. ( 49 ) 

Because g 2 = and h 2 = h 1 we are not compelled to use all of Equations 40 to 45, but only 
four of them, say, Equations 40, 42, 43, and 45. 

We consider the main problem in the theory of satellites: to determine the motion of the 
satellite under the influence of the disturbing function 


m 3 r 2 

— — p 2 (*) 

r'3 


( 50 ) 


where a is the cosine of the angle between r and r'. In the beginning we omit the parallactic terms 
of the first and higher orders. Their inclusion with accuracy up to k 3 does not present any 
difficulty. It is necessary to emphasize that in any attempt to achieve a higher degree of accuracy 
it is meaningless to include the higher harmonics alone without taking the "cross actions" between 
the lower harmonics into consideration. 


Taking into account Kepler's law, 


nij + m 2 + m 3 
m i ^ m 2 

we can write Equation 46 in the form 




n = 


m 


3 


m 2 m 3 






P 2 (o') dl. 


(51) 


If e # not small, then in treating the long-period effects it is better to use f ', the true anomaly of the 
distant companion, as the independent variable instead of time. This idea was suggested by Brown 
(Reference 3) in his work on the stellar three-body problem. 


In accordance with the form of the differential equation of motion of the satellite given in Equa- 
tion 5, we have 


or 


df' 

~dt 


Oj + m 2 + m 3 ) a 1 (l - e ; 2 ) 


.i 2 df' rV^-3/2 
dt n 


r j 2 a * 3/2 (l-e' 2 ) 


- 1/2 2 


(52) 


Putting 


5 = exp i f' 
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and introducing the operator 


D = l 


d £ d(if') 


(53) 


analogous to the operator 


d(n - n')i t 

of Hill's lunar theory, the system of Equations 40 to 45 becomes 


^ „ / BF , BF \ /, BF BF 

Dgl = 2 ( g3 air + 3 — ) + ( h i — - 




3g. 


3h, 


(54) 


3F , 3F 

D e 2 = 2 g 3 ^r- - h 
Bh, 


— ) * K 5L 


+ g~ 


BF 

BhI 


(55) 


BF 


BF 


Bh 


Bh 


D «3 = ( g i +g 2 ^r)- h i BTT" +h 2 


BF 


1 


+ h 


BF 




(56) 


Dh, = 2 (h 


(- s, ■ a - 1 


l , 3F BF 

“ I h + g, 

1 ^3 3g 3 


(57) 


and 


where 


and 


Dh a = -2 h 3 


BF BF \ / 

+ g 3 ] + ( h 

3h i 3 V 


3F . jF 
2 2 


,3 F 3F \ /, 3F . BF 

Dh 3= ^--*a — W h i — " h 2 — 


Bgj 


Bh 


Bh, 


F = \(1 + e' cos f 


■> ■ d ( (!)■ p 2 <-> 


3 ^ (1 - e' 2 )"3/2 


m ! + m 2 + m 3 n 


(58) 


(59) 


(60) 


(61) 


is Brown's parameter. The use of k permits one to write the final results in a more compact form, 
thus avoiding a considerable portion of the development into the power series relative to e'. 


In his earlier work (Reference 12) the author obtained the formula 


f 2?r 

L. J Q 2 P 2 (<y) di = i [ 5(e • u 0 ) 2 + (C X n 0 ) 2 - e 2 1 


(62) 
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where u° is the unit vector in the direction of 7' , and e and c are defined by means of Equations 3 
and 4*. We have 


and 


e • u° = cos f' + e 2 sin f' = -i-(gj £*‘ + g 2 £)> 

c • u° =i(h 2 c-hj r 1 ). 


(63) 


(64) 


Substituting these values into Equation 62 we obtain, after some easy vectorial transformations, 
taking Equations 47 and 48 into account; 

F=-^-\(l + e ' cos f<) [ (5g i 2+h i 2) ^ 2 (65) 

+ (6 g, g 2 + 2 hj h 2 + 4 g 3 2 + 4 h 2 ) 

+ (5 g 2 2 + h 2 ) £ 2 ]. 

Substituting this value of F into Equations 54 to 59 we obtain 

Dg, = J Mi + e' cos f') (a a 0 + a 1+2 £ 2 ), (66) 


and 


where we put 


0g 2 = -~K(l + e* cos {') (/3, _ 2 r 2 + /3 li0 ), 


(67) 


Dg 3 = — \(1 + e' cos f') ( 7 j _ 2 r 2 +y li0 + r 1+2 5 2 ). 


( 68 ) 


, Dh l = f ' K(1 + e COS f ’> ^ a 2 , 0 + a 2 , +2 t )' 


(69) 


Dh 2 = -~\(1 + e' cos f') (fi 2 -2 Z + ^2,o>- 


(70) 


Dh 3 = g-k(l + e' cos f') (y 2 _ 2 r 2 + r 2i+2 C 2 ). 


(71) 


a i,o “ ^ g 3*h + 


a i, +2 = 5 & 2 h 3 + g 3 h r 


(72) 

(73) 
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Pl,-2 ~ 3 E 1^3 “ g 3*V 

(74) 

A.O = - 3 g 3 h 2 + g 2 h 3> 

(75) 

>1,-2 = - 2g l h l* 

(76) 

> 1.0 = - g i h 2 - 

(77) 

>1 , +2 = ~ 2& 2 h 2 , 

(78) 

a 2 , 0 = ~ 3g l g 3 " ^1 ^3' 

CD 

a 2, + 2 = - Sg 2 g 3 + h 2 h 3- 

(80) 

^2.-2 = 5g l g 3 + h l h 3* 

(81) 

02.0 = 5g 2 S 3 _h 2 h 3’ 

(82) 

>2.-2 = 5g ? + h J- 

(83) 

>2 , +2 = - Se 2 - h 2' 

(84) 


KRYLOV-BOGOLUBOV METHOD 


The Krylov-Bogolubov method (Reference 7) has not yet achieved general recognition in lunar 
problems: A recent application of this method by Morrison (Reference 13) to the planar lunar 
problem deserves especially to be mentioned. We shall give here a short outline of the method of 
Krylov-Bogolubov from the viewpoint of its application to the stellar three-body problem. The 
differential equations we are considering here are of the form; 


Dgj = \Gj (gj, g 2> g 3> hj, h 2 , h 3 ; £; \), 

(85) 

Dg 2 = ^ G 2 ( g l> g 2’ g 3’ V h 2’ h 3 : 

(86) 

Dgj = (g t > g 2 > g 3 > hj, h 2 , h 3 , £; \), 

(87) 

Dhj = \Hj (g (> g 2 , g 3 , hj, h 2 , h 3 ; £; K), 

(88) 

Dh 2 = XH 2 (gj, g 2 , g 3> h 2 . h 3 ; £s X), 

(89) 

Dh 3 = ^ H 3 g 2 . g 3 . h x , h 2 , h 3 ; £; X), 

(90) 


where g i? G 2 , G 3 > H, , H 2 , and H 3 are the series of polynomials in £ and 1 with the coefficients 
being polynomials in g„ g 2 , g 3 , h p h 2 , and h 3 . 
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The Krylov- Bogolubov method consists in finding a transformation 

g l = g l + (gj, gj, gj, hj, hj, hj; £; A), 

g 2 = g* + Abj (gj, gj, gj, hj, hj, hj; £ ; A), 

g 3 = g* + (gj, gj, gj, hj, hj, hj; £; A), 

hj = hj + Aa 2 (gj, gj, gj, hj, hj, hj; £; X), 

h 2 = hj +X.b 2 (gj, gj, gj, hj, hj, hj; £; X), 

h 3 = h 3 + ^ C 2 < g J> g 2’ 63’ H ?’ h 2’ h 3 ’ ^ 

so that the differential equations in the new variables do not contain £ : 


and 


and 


Let us introduce the operators 


K = £ — 

3£ 


and 


A = GJ — + g; — + Gj — a 


+ Hj — + H 2 — + H; 


3g 2 ° 3gj 1 3h l on 2 “‘‘3 


d 

3h! 


d_ 

3 3hj 


'1 , T ~2 

3gj 

Then the operator D as applied to Equations 91 to 96 can be written in the form 


D =K + XA. 


It follows from Equations 91 to 105 that 

Dgj = AGj + A.Kaj + X 2 A a^, 


(91) 

(92) 

(93) 

(94) 

(95) 

(96) 


Dgj = AGj (gj, gj, gj, hj, hj, 

1*3’ 

(97) 

Dgj = AG* (gj, gj, gj, hj, hj, 

h- A), 

(98) 

Dgj = AG* (gj, gj, gj, hj, hj, 

hj; A), 

(99) 

Dh* = AH* (gj, gj, gj, hj, h* 2> 

hj; A), 

(100) 

Dh* = AH* (gj, gj, gj, hj, hj, 

hj; A), 

(101) 

Dhj = ^ H 3 §3- hj, hj, 

hj; A). 

(102) 


(103) 


(104) 


(105) 


(106) 
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and 

Introducing the shift operator 
1 + T = exp 


Dg 2 = + \Kbj + A 2 Abj, 

Dg 3 = AG* + \KCj + A 2 Acj, 
Dh, = AHJ + \Ka 2 + A 2 Aa 2 , 
Dh 2 = AH* + \Kb 2 + A 2 Ab 2 , 


Dh 3 = AH* + AK c 2 + A 2 A c r 


( a . 3 

I a i + b i 

\ * ^gj 3 g 2 


s 2 ^g; 

and combining Equations 85 to 90 and 106 to 111, we obtain 


3 3,3 3 

+ c, + a + b, + c 

3 ^ 3h 2 ^*3 


and 


(107) 

(108) 

(109) 

( 110 ) 
( 111 ) 

( 112 ) 


a; + AKaj + A 2 Aaj = d + T) AG r 

(113) 

^ +AKb t +A 2 Abj = (1 +T) AG 2 , 

(114) 

3; + AKcj + A 2 A Cj = (1 + T) AG 3 , 

(115) 

a; + AKa 2 + A 2 Aa 2 = (1 + T) AH,. 

(116) 

«* + AKb 2 + A 2 A b 2 = (1 + T) AH 2 , 

(117) 

h; + AKc 2 + A 2 Ac 2 = (1 + T) ah 3 , 

(118) 

3 g,, g 2 , • • • are replaced by g * , g* 2 , . . . . 

, bj , . . . c 2 , and G* , G* , . . . H* in the form of power series in , 

\a 1 = \ a n + k 2 a 12 + . . . , 

(119) 

\b 1= \b 11+ x* bl2 + ... , 

(120) 

^ c i = \ c n + \ 2 c i 2 + . . . t 

(121) 

\a 2 = \a 21 + k 2 a 22 + . . . # 

(122) 

A. b 2 ~ A.b 2 ^ ^ ^22 + * * * > 

(123) 

A- C 2 ~ ^ C 2l + ^ C 22 + * • * » 

(124) 
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\G* = A-G'j + A 2 G * 2 + . . . , 

A.G* = A.G*j + A 2 G * 2 + . . 

^g; =w ? 1 +x 2 g * 2 + . . . , 

AH* = \H*j + X. 2 H * 2 + . . . , 

AH; = AH 21 + A 2 H * 2 +. . 

and 

**$ = **; i +^ h ; 2 + — 

We have similar developments for the operators 1 + t and AA, 

1 + T = 1 + ATj +\ 2 T 2 +. . . , 

the operators T k are Faa de Bruno differential operators (Reference 14) 

AA = AAj + A 2 A 2 + . . . , 

where, from Equations 112 and 104, we deduce 


T 


l 





+ 



(125) 

(126) 

(127) 

(128) 

(129) 

(130) 

(131) 

(132) 

(133) 


A 


l 



+ G n 



H u 





(134) 


From Equations 113 to 130 we obtain 


G 32 + K C 1 2 + ^1 C 1 1 - T 1 


'll =<v 

H n 

+ Ka 21 

= H r 

cs 

a 

ii 

H 21 

+ Kb 21 

= h 2 - 

ii 

o 

Co 


+ Kc 2 i 

= h 3 . 

= T 1 G l- 

h; 2 

+ Ka 22 

+ Aj a. 

II 

k- h 

o 

to 

h ; 2 

+ Kb 22 

+ Aj b. 

- T G , 

1 3’ 

h ; 2 

+ Kc 22 

+ \ c. 


(135) 




(136) 
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•Civ 


Let us make use of the Krylov- Bogolubov averaging operator M. In our case its role consists 
in removing the positive and the negative powers of £ from the expression to which it is applied, 
thus reducing this expression to the term which is independent of 5. In addition to M, it is convenient 
to introduce the operator P by extracting the nonzero powers of £ from the given expression. It is 
useful to remark that 

MK = 0, 


and 


M P = PM = 0, 


MK" 1 P ^ 0, 


(137) 

(138) 


where K -1 designates the inverse operator of K, and the operators A, A,, A a , . . . , and M are 
commutative: 


MA = AM, 


MA. = A. M, j = 1, 2, 3, • • • 


(139) 


From Equations 138 and 139 we have, in addition, 

MA IT 1 P = 0. 

Taking Equation 137 into consideration we deduce from Equation 135: 


G » = MC V 


G » = MG 2- 


G « = MG 3' 


Hj, = MHj, 


h* 2 , =mh 2 . 


h ' 31 = mh 3 , 


(140) 


(141) 


a u =K-> P Gj> 

b u =r» pg 2 , 

C U = rl PG 3- 

From Equations 137, 138, 140, and 142 we obtain 

G i 2 = MT i G r 
G 2 2 = MTj G 2 , 

G «= MT . G 3> 


a 21 = K” 1 PH,, 


b 2 , = K-i PH 2 , 


C 2X =K->PH 3 . 


H* 2 = MT, H,, 
H* 2 = MT, H 2 , 

H 32 =MT , H 3* 


(142) 


(143) 
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and, as a consequence, we have 


a 12 = K -1 PCTjGj - Aj Sjj). 

b i2 = K_lp (T 1 G 2 -A x b u ), 
c 12 = K-i P(T 1 G 3 _A 1 c u ), 


a 22 =K-i P(T H x - Aj a 21 ), 

b 22 =K-i P (Tj H 2 - Aj b n ), 

C 22 = ^ ^^1 H 3 “ ^1 C 21^’ 




> 




(144) 


This process is continued until all the significant terms containing £ are eliminated. 

The differential equations for the elements affected only by the long-period terms become 

D g ; = M(1 + T : ) Gj + . . - , Dh; = M(1 + Tj) Hj + . - . , ^ 

Dg; = M(1 + T % ) G 2 + . . . , Dh; = M(1 + T x ) h 2 + . . . . > (145) 


Dg; = M(1 +TP c 3 + . . . 


Dh; =M (1 + T x ) H 3 + . . . . 


J 


In further exposition we shall change the notations and shall write g 1? g 2 , . . instead of gj, gj, • . • . 

INTEGRATION OF DIFFERENTIAL EQUATIONS OF STELLAR THREE-BODY PROBLEM 


From the Equations 66 through 71 we have: 


G 1 if (1 + e ' cos n ^ a i,o + a i,+ 2 £ 


+ 2>> 




G 


2 


--J(l +e' cos f‘) (p u _ 2 T 2 +/ 8 lf0 ), 


G 3 = i(l + e' cos f') (y u _ 2 C 2 + 7 , 


1 ,0 


+ y li+2 £ +2 ), 


> (146) 


--(1 + e' cos f') (a 2 0 + a 


2 , +2 


£ +2 ). 


H 2— jd 


+ e' cos f') 03 2 '_ 2 r 2 


+ ^2,<P’ 
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H 3 = -|(l + e' cos f'l (y 2 . 2 T 2 + 7 2 , +2 £ +2 ). 



Making use of Equations 141, 142, and 146 we obtain: 


G n =+T a ,,o' 


G * 21 =~zAo’ 


Tf* ° 

H 11 " + 4 a 2,o 


H* 

21 


■k. 


G ;i=+T^.O- H 31 =0 • 


(147) 


and 


a U = |" a i ,+2 ^ + "I e< ^1 ,0 T 0 + a i,+2 T +2^ 


b ll = t/^1,-2 ^ 2 -“ q " 6 ^ 1,-2 T -2 + ^ 1,0 T 0 ^’ 


= -|-(7 lt _ 2 r 2 -7 1>+2 ^ 2 ) + fe' (7 1( - 2 t - 2 +yi.o T o+^. + 2^ 


3 21 = !T a 2,+2 ^ + g" e ' ( a 2,0 T 0 + a 2,+2 T +2^ ’ 


r ( 148 ) 


b 21 “ o ^ 2,-2 ^ 2 o e ^ 2,-2 T “2 + ^ 2,0 T 0 ^’ 


C 21 “ “ 16 ^2,-2 ^ 2 ^2, +2 ^ ^ 


— e' (7 2 - 2 T “2 + ^2. +2 T +2^’ 
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where we put 


t -2 = - J r3 - r1 ’ 


T o =-t _I +C- 


> (149) 


T *7 = + 5 +4- £ 3 
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From Equations 133 and 145 to 148 and taking 


M (£~ 2 t - 2 cos f') = 0 » 

Mcr 2 cos f') = 

M(r 2 r +2 cosf')=|-, 
into consideration, we obtain: 


M (r 2 cos f') =-y> 
M(r Q cos f') = 0, 
M(t + 2 cosf') = y, 


M(£ 2 r cos f') = -— , 
3 

M(£ 2 r Q cos = 

M(£ 2 t +2 cos f 1 ) = 0 


Dg i = | Xa l,o + |^ 2 a +|e* 2 ) (- 5g 2 7 2> . 2 + 2g 3 /? 2i _ 2 - 2 h 2 7 j, _ 2 + lOh^ _ 2 ) 


, -^ 2 e' 2 
+ 64 


+ ^ gj (>2, +2 “^2,-2^ + g 3 ^2,0 + 3a 2, + 2^ 


+ 3hj + 2 -7 lt _ 2 ) “^2 ^ 1,0 + ^3 ( a i, +2 + ^^i f o) J ’ 


Dg 2 ~ 4^1,0“ 64 ^ 2 ( 1 + 3 e 2 )( +5g 1^2.+2 2g 3 a 2, +2 “ 2h l'> / l,+2 + 10h 3 a i,+2) 

-fr* e#2 [ + f g 2( 7 2,+2“ r 2,-2)" g 3( a 2,0 +3 ^2,-2^ 

+ 3 h 2 (7^ _ 2 ” 7j ? + 2 ) “ h ! ^1,0 + ^3 , - 2 +5a lt o)] » 

Dg 3 = f x r ll0 - ^K 1 + 3 e ' 2 )< g i a 2, +2 + g 2 ^ 2> _ 2 + hja^ +2 + h 2 p u _ 2 > 

“(54 e [ g l(2 a 2l o +^2, - 2 ) +g 2 (^^2,0 +a 2, + 2 ) 

+ hj (2“ 10 +<®j, - 2 ) + ^*2 (^^ 1,0 + °i, + 2 ^ J ’ 


Dh i =1^2.0+^ ^ 2 ( 1+ 1 e ' 2 )< +10g 2^.-2-10g 3 /8 lt . 2 -h 2 r 2i _ 2+ 2h 3 ^ 2> _ 2 ) 
+ ^- 2 e 2 £ + 5gj CXj-2 -y li+2 ) + 5g 2 7 10 -5g 3 (^1,0 + a i, +2^ 

+ "2 ^1 O 2 , '2 - ^2, + 2 ) + ^3 (^2, 0 “ a 2. + 2 ) 1 ’ 


(150) 


(151) 


(152) 


(153) 
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Dh 2=-f^ 2 .-2 


" e ' 2 )<+ 10g l^l. +2 + 10g 3 a i. +2 +h l> , 2. +2 + 2H 3 a 2. + 2> 


e ' 2 [ + 5g 2 Cy 1>+2 -r 1> _ 2 ) + 5 g 1 ri,o + 5g 3 (a 1 ,o+^i,-2) 


+ ~2^ 1 2^ , 2.~2 ~^2, + 2 ) + * 1 3( a 2,0 “^2,-2)] ’ 

(154) 

Dh 3 = ^ e ' 2 [ 5 (e, a i,0 * 6 2 / 3 1 ,o) + ] • 

(155) 

h 2 = h , we are not obliged to use all of Equations 150 to 155. 

The choice of 


four, Equations 150, 152, 153, and 155, will suffice. 
The formulas 15 to 20 suggest the substitution 


gi = e Pr g 2 = e P 2 > g 3 = 6k P 3 ’ 

h j — k w j, h 2 =kw 2 , h 3 =w 3 , 

Here p x , p 2 , p 3 , w 15 w 2 , and w 3 are the functions of £ Cl , C 2 and e 2 , k 2 ; c t is the mean motion of 
the longitude of the perigee; and c 2 is the mean motion of the longitude of the ascending node. For 
the values of e and k small enough, the quantities g x , g 2 , g 3 , h 1? h 2 , and h 3 are developable into 
series in eC l , e£ Cl , k£° 2 and k£ ° 2 ; both cj and c 2 are series in e 2 , k 2 , e' 2 , and a 2 , where a. is the 
parallactic factor. We consider here the case of a close satellite and do not discuss the influence 
of the parallactic terms. The theory presented here could be extended further by including the 
higher powers of k and a. 

The fact that the Moon is a body and not a point, together with the assumption that the paral- 
lactic factor is small, imposes a limitation on the upper value of e. We shall consider the case 
when k is not small, but let us assume that e is small or moderate, say e £ 0.3, or less. If we 
make use of the substitution (Equation 156) the differential equations become: 

Dw l =~f Xw i w 3 e' 2 )(2w lW |-k2w2w 2 ) 

+ \e 2 A 1 +i 2 e' 2 B 1 , (157) 
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(158) 


lllll■llllln 


Dw 3 =\ 2 e' 2 B 3 , 

Dp l = f’*-(Pl w 3 + 3k2 p 3 w l) 

+ ^-\ 2 (l +| e' 2 ) [50p 1 w 2 +k 2 (-5p 2 w 2 -8p 3 w 1 w 3 +4p 1 w 1 w 2 ) ] 
+ \e 2 F 1 +A. e' 2 Gj , 


and 


DP 3 =-f^ (P l W 2 + P 2 W l) 

-21 X 2 (l+| e ' 2 )w 3 (pjW 2 tPjWj) +\ e' 2 Q 3 , 


where we put 


. IS 
A ! = -T p i p 3 


+ ^A(l + f e' 2 ) (-20p 1 p 2 w 1 -5p 2 w 2 -40p 1 p 3 w3^+10k 2 p 2 w 1 


B -- 2 - 

1 64 


-2w 2 w 2 + |k 2 w^w 2 +w 2 J 


+ -L e 2 

64 


-4 P ? W l-| P 2 W l +5P l P 2 W 2- 2 ° P 2 P 3 W 3) +10k2p l W 2 


B 3 = S k SH' w '0 + i £2w 3( p ?- p ^) 


+ 32 e2k2(PlWl+P2W2> ’ 


F i = |f( 1+ 5 e' 2 )(-25p 2 P 2 + 10k 2 p lP : 


(159) 

(160) 

(161) 

(162) 

(163) 

(164) 
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G 1 - 



w 


2 

3 


and 


+ — k ^( + U pi V 


• 24P3W2W3-10PJWJW2 +pjW2), 


(165) 


27 


Q 3 = --^ W 3 ( r >lW l + P2 W 2) 


+ ~e 2 P 3 <P\ -P|) -~k 2 P 3 (wf -«!). 


(166) 


In the lunar theories of Delaunay, Hill, and Brown the motion of the node and of the perigee 
are developed into the power series in e 2 , k 2 , e ' 2 and a 2 . However, the domain of covergence 
of these series is not given and their application to the cases of large inclinations and large ec- 
centricities cannot be taken for granted. The choice of the small parameter k and the appearance 
of the closed form factor 1 + (2/3) e' 2 remove the necessity of development relative to e ' , as 
it has been found by Brown. In order to remedy the situation relative to e 2 and k 2 , we suggest 
here the numerical computation of c t and c 2 by successive approximation rather than development 
into power series. The numerical values of e, k, e' , and X are substituted from the outset. 

We have to rewrite Equations 157 to 160 and also to develop the equations for the determina- 
tion of cj and c 2 in a form suitable for the application of the process of successive approximations. 
As a by-product we shall also establish the domain of convergence to the inclination. 

The integrals 34 and 49 take the form 


p i w 2 ~ p 2 w i + 2 p 3 w 3 = 0, 

^k 2 WjW 2 + e 2 ^p 1 P 2 - k 2 p 2 ^ = 1 . (168) 

The last equation is to be used for determining x 0 at each cycle of the iterative process. The 
first equation can serve as a check of computations. 

We have three groups of terms in the right side of Equations 157 to 160: the terms propor- 
tional to those proportional to k 2 (1 + (2/3) e' 2 ), and those proportional to k 2 e' 2 . The last 
two groups were produced by the mutual action of the short-period terms. The terms of the first 
two groups correspond to the terms included by Brown in the theory of the stellar three -body 
problem. They are the source of secular changes in the argument of the perigee and the longitude 
of the node, and in addition they produce the main long-period effects in the elements. In the per- 
turbations of canonical elements of Delaunay, the first two groups give rise to the terms depending 
upon the mean argument of the perigee. The terms proportional to X. 2 e /2 are omitted in Brown's 
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theory. To the accuracy with which we are concerned, they have no influence on the mean motions 
of the arguments; they will produce in e and c only the long-period effects proportional to \ e ' 2 . 


We have 


and 


where we set 


DWj = Xv n w t + A. W r 


(157’) 

Dw 3 = K 2 e' 2 W 3 , 


(158’) 

D(£ C2 P x ) = V n (£ C2 P 1 )+^ 2 M 12 k 2 (£ +C2 p 2 )+\ Ml3 k 2 P 3 + ° 2 

Pi)- 

(159’) 

Dp 3 = ^31 [(£ 1=2 Pi) + (£ + ' 2 Pj)J + ^ P 3’ 


(160’) 

v » = -f - k ’>' 


(161’) 



(162’) 

--If (’*!«'•)■ 


(163’) 

9 9 . /, 2 ,9 ^ 3 27 . ( - 2 , 2 

^3= IT “ ¥ v 1 + 7 e / x o’ ^31 = - T - iH 1 + T e 

) *0 ’ 

(164’) 

Wj = e 2 Ai - A. (l +| e’ 2 j k 2 w x (w, w 2 - 1 ) 



"f W l (W 3 “ *c) +Xe ' 2 B 1 ’ 


(165’) 

w 3 =b 3 , 


(166’) 

P x = +-- k 2 (w x - C° 2 ) P 3 



" 64 ^ (* + T e ’ 2 ) k2 [ 5(W2 " S+2C2) P2 + 8 X o <>i ~ ^ +C2 ) p 3 

] + e 2 Fi 


+ -J p i ( w 3 " V + ^ e ' 2 G x . 


(167’) 
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and 

P 3 = M 31 [(w 2 - C’ C2 ) Pi + (w t - C 2 ) P 2 ]+ A. e' 2 Q 3 . (168’) 

Eliminating p 3 from Equation 159’ by means of Equation 160’, we obtain 
D 2 (£' Ca Pj) = \ D(£ Ca Pj) + >^ 2 m 12 k 2 D ($ C2 p 2 ) 

^ 2 M 13 M 31 k 2 [a' 2 p,) + (C +C2 p 2 )] + ^ Sl , (159”) 

where we set 

\ 2 S, = X 2 Ml3 k 2 P 3 + X D (£‘ C2 Pj) . (159™) 

Equations 157*, 158 r , 159", and 160 T are in a form suitable for the application of the method of 
successive approximations. 

In order to describe the process of integration of Equation 159 T ’ let us consider a typical 
term and a typical differential equation 

D 2 (i ° 2 Pj) = K m u D(£ ° 2 p t ) + A 2 m 12 k 2 D(^ 2 p 2 ) 

+ X 2 / x 1 3 M 3 1 k 2 [(c' C2 p 1 ) + (C° 2 P 2 )] + ^ 2 (NjC a +N 2 r a ) . (169) 

where a = (2k 1 + i) c 1 + (2 k 2 - 1) c 2 . 


Using the substitution 


T 2 Pl = k x c a + k 2 c a , 


we obtain 


or 


Setting 


and 


a 2 K, = \ a K t + \ 2 m 12 k 2 a K 2 + X 2 m 13 M 31 k 2 (Kj + K,) + X 2 N t = 0 , 
a 2 K 2 = - K Mn a,-X 2 m 12 k 2 a Kj + X 2 m 13 M 3l k 2 (K, + K 2 ) + \ 2 N 2 = 0 , 

(a 2 - 2 X 2 k 2 Mi3 /x 31 ) (K, +K 2 ) - a <\/x n - \ 2 k 2 Ml2 ) (Kj -K 2 ) = X 2 (N, + N,) , 
a(\ Mn + \ 2 k 2 m 12 ) (Kj + K 2 ) - a 2 (Kj - K,) = ^ 2 (N 2 - Nj). 

c l = k er x . 

c 2 = A. cr 2 , 
a, - Cr 2 = a , 

/8= (2kj + 1) a, + (2 h 2 - 1) o-j, 


(170) 
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where cr x and a 2 are of the zero order relative to e, we deduce from Equation 170 
(/ 3 2 -2 k 2 Ml3 m 31 ) (Kj + K 2 ) - p (M u - Xk 2 Ml2 ) (Kj - K 2 ) = N 2 + N t , 


/3(m u + \ k 2 Ml2 ) (Kj + K 2 ) - /3 2 (Kj - K 2 ) = N 2 - Nj . 

The determinant of this system is 

A = P 2 (Ml! + 2 k 2 Ml3 m 31 - X 2 k 4 m 2 2 - /3 2 ), 
and the determinants associated with + k 2 and - K 2 are 

- /?(M U - ^k 2 M 12 + P) Nj + /3 (/z u - kk 2 fx l2 - N 2 , 
A 2 = - /3(/x n + \ k 2 /i 12 + p) Nj - /3 (m u + kk 2 - /3) N 2 


(171) 


(172) 


(173) 

(174) 


+ 2 k 2 ^3 M 31 (Nj N 2 ). 

The solving of this system is a straightforward process and, generally speaking, the division by A 
does not introduce any small divisor. 

For 

*1 - *2 = “ 1 » 


the part independent from 6 and k disappears in j3 . For small values of e and k the value of c l + c 2 
then will also become small. Thus, this case requires a special discussion, and we shall show 
that in the problems of integration, even for small values of 6 and k, no small divisors will ap- 
pear which would make the very coefficients of £ j(c i +c 2 ) anc j ^ j < c i +c 2 > i ar g e ^ 

A special consideration is also required for the case 


because it is connected with the determination of the mean motion of the pericenter and with the 
finding of the first approximation to . 

In agreement with the theory of the long period effects by Poincare (Reference 4), we have the 
following form of solution of £ 2 Pl , P 3 and of £~ 2 Pl , p 3 : 


,(2i +l)c + (2i 2 -l)c 

1 1 z z + A e 


' 2 Y! Mi v i 2( £2 < k2 ) £ <2i i +1)c » + (2j 2' ° c 2 ; 


£ N illi2 (* 2 . k 2 H <2il + 1)C ‘- 
and of £" c 2 Wi andw 1 : 

K i,.i 2 ( e2 ' k2 > S 2ilCl + 2l2C2 + Xe ’ 2 ^] L i,.i 2 ( €2 ' k2 ) ( 2 * 1 C ! + 2 ^ 2 C 2 ; 
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and of 


3 - * 0 


In addition 


- x and X.e # 2 w 3 : 


and for p 3 , 


Xe ' 2 Z] Li i’ i 2 (e2> k2 K 2il ° 1 + 2i2C2 ' 


111 + i 2 l = ° . 


1 3 1 + J 2 I = 1 ; 


N. . = - N_ . +1 

1 l* l 2 i!-!. 1 2 +1 


M. . = -M_. . , ; 

1 1 ’ 1 2 L 1 l ' 1 2 + 1 


but for P a 


N. 

i 

M. . 


, , i~ — N . . 

12 1,1, i 2 +l 


^•‘2 - M -i,-i,-i,+r 


n x 2 


In our case we have to consider the terms 


C 1+Cl and £ 


: i c 2 


Retaining only these terms in P 3 and £ C<1 P 1? we can set 


and 


P 3 =M : £ Cj C2 + M 2 f Cl " 


C 2 ^ =Cj C * 1 '* 2 +C 2 C l ~ 


2 


2 


where M 1? m 2 , C 1? and C 2 are the functions of e 2 and k 2 . If e and k both are small, then we shall 
assume that they are of the same order of magnitude. 

We have to set 


N, = M 13 k 2 M, + /3C, , 
N 2 =M 13 k 2 M 2 -0C 2 . 


Then Equations 172 to 174 become 

A = ^ 1 /3 2 +. . . , (175) 

Aj = -/3 2 Mll (C, +C 2 ) +k 2 /3 Mll Ml3 (M 2 _M,) + . . . . (176) 
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A 2 = - (C t - C 2 ) + 2k 2 /3^3 m 31 (C t + c 2 ) 


(177) 


- 2k 4 ^ ^ (M, - M 2 ) - k 2 /3^ 3 (M, + M 2 ) + . . 

where only the terms of the lowest order are written explicitly. We have 

K, + K = - — (C. +C 2 ) +^J— (M 2 -M ) + • . . , 
/'ll P /'ll 


(178) 


v V 1 /p p \ k /'l 3 /^31 /p p» \ 

K 1 - K 2 = -— (C 1 + V + ~B ^ (C 1 + ^ 

M X1 ^ Ml! 


(179) 


k 4 , 2 /'l 3 ^31 


11 


(M, _ M 2 ) -jj. (M, + M 2 ) + 

P Mi t 


From differential Equations 157 to 160 we easily deduce as the first approximation for small values 
of k and e : 




-M* 4 - ! ) 


+ c 2 = - — \(3k 2 + 5 e 2 ) + . . . , 


Consequently, 


n » , 2 13 2 

fl = k' 

H 8 8 


kf. 

/3 


and from Equations 178 and 179 we can see that in the process of integration the order of the coef- 
ficients of long- period terms satisfying the condition i a - i 2 = -1 remains the same as that before 
the integration. 

In the case a = c - c 2 , we have 

K i = I* 

P = °i - °2 = a 

and Equations 170 become 

(O- 2 - 2k 2 Mll ) (1 + K 2 )-a (Mn - M< 2 Mj2 ) (1 - K 2 ) = N 2 + N,. (180) 

cr( Mu + k 2 Ml2 ) (1 + K 2 ) - a 2 (1 _ K 2 ) = N 2 - Nj . (181) 
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In order to simplify the elimination we add the expression 


Cl +K 2 ) +(1 -K 2 ) =2 . (182) 

The coefficients N 1 and N 2 are taken from the corresponding part in s x . 

Elimination of K 2 leads to 


^ - 2k 2 M 12 ^31 

-*<*11 - k2 ^12 ) 

n 2+ n, 



+ k 2 M 12 ) 

-0 2 

N 2- N l 

= 0 

(183) 

1 

1 

2 




for the determination of o- and, consequently, of c 1% Equation (183) can be written in a form 
suitable for the application of the method of successive approximations: 

0-2 “ + 2k2 ^13 M 31 - k4 + r ’ (184) 

where we put 

^ k2 Mi 2 \ 

1 J (N 2 + n x ) 

+ (N, -N 2 ) ( Mu - k2 M 12 ^31 )• (185) 

The terms contained in r are of the higher orders. They contain either e 2 or k as a factor. The 
value of cr in Equation 185 can be taken from the preceding cycle of iteration. 

To obtain the first approximation to c x and c 2 in the case of a large inclination, we neglect 
e 2 and k 2 and set 


3_ 
4 


^ - k2 • 

In accordance with Equations 162 f to 164* and 184 we have 

a i - tr 2=T /-| k2 • 

Thus 


(186) 


(187) 


(188) 


27 



(189) 


C 2 = -|\ x/l - k 2 . 


We see that for a small value of e we have a resonance effect at k 2 = 2/5, I % 39°, or I % 141°. The 
critical argument is the argument of the perigee. 


We have to determine the value of K 2 . From Equation 181 we obtain 

- ^ k2 M 12 + [( N 2 " N i^ /cr ] 


K, 


oh + M u + ^ k2 M 12 


In the first approximation we have 

a -^u 3 


K 2- 


k 2 


<r + M, 


(A 7 !* 7 * ^ Tir ’) i 


and 


P, = S +Cl +k 2 f Cl+2 ° 2 , 
p 2 = f Cl + K 2 C*' 1 "*' 2 , 

Wj = i +C2 , 


w 2 = £ t2 , 


(190) 


(191) 


(192) 


where K is defined by means of Equation 191. Substituting these values in Equation 160’, we have 


Dp , = (1 + K ) ( £ 1 ° 2 + £, 1 2 ) . 


and taking Equation 187 into account we obtain the first approximations to p 3 : 
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Using the solutions in Equations 192 and 193, we start the process of iteration and continue it till 
we reach the convergence. The integration of Equations 157’ and 158' is a straightforward process 
and requires no additional explanation. 


COMPARISON WITH DEVELOPMENTS OF DELAUNAY 

The lunar theory of Delaunay gives the mean motions of the argument of the perigee g and 
of the longitude of the ascending node h in the form of power series in the mean values of n'/n, 
y - sin 1/2, e, e' and a' /a. The information about the domain of convergence of these series 
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cannot be obtained in an easy way. We shall attempt to obtain some information about the con- 
vergence using the closed-form expressions for ( Cl - c 2 ) 2 and c 2# Delaunay*s development of 
dg/dt corresponds to c x - c 2 of our theory; dh/dt corresponds to c 2 . We consider those parts of 
dg/dt, dh/dt, ,c 1? and c 2 which do not contain e because previously we made the assumption that it 
is small. We mentioned also that if this assumption is made then the solution can be obtained in 
the form of trigonometric series with the usual arguments of the lunar theory. 

From Equations 157* and 184 we have accurately up to A. 3 , 

( c i - c 2 ) 2 = X 2 (^\ x + 2k 2 /x 13 ^ 31 )> (194) 

and, accurately up to \ 2 , 

c 2 (195) 


If £ = 0 then we have k 2 = 1 - x 0 2 , and making use of Equations 161* to 164* we obtain with 
the same accuracy as before the closed-form expressions 


< C 1 - C 2^ 2 =| ^ 


and 


.'>)( 25 x ;- 9 ) Xo ' 


(196) 
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The results 196 and 197 can be extended to include higher powers of A and of the parallax. Of 
course, the parallactic factor must be properly modified. The solving of this problem will consti- 
tute the topic of future work. It will require keeping the higher orders of the **cross actions’* of 
short- period terms and parallactic terms. In order to facilitate comparison of our results with 
those of Delaunay, we shall eliminate x 0 from Equations 196 and 197 by means of 


*0 = 1 - 2 72 • 


and perform the expansion in powers of y we obtain 


< c i - c 2> 2 = f 


and 


j(l - 10 y 2 + 10 /) + i A. (l+| e' (4 - 33 y 2 + 75 / _ 50 y 6 )| 


c 2 = ~^k(l - 2-y 2 ) + -^X. 2 + § e * (1 - 6 y 2 + 6 y 4 ) . 


Assuming that y is small, we obtain from Equation 198: 


c i - C 2 = 4 




(198) 


(199) 


( 200 ) 


keeping only the terms necessary to make the comparison with the results of Delaunay. 
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In the definition of A. the factor (l - e' 2 )~ 3/2 is included. The presence of this factor in the 
expanded form in Delaunay's results is beyond any doubt: it can be seen directly if one looks care- 
fully at the developments of dg/dt anddh/dt. 

Thus, we have to compare the parts multiplied by \ and \ 2 with the corresponding y -terms 
multiplied byn'/n and (n') 2 /n 2 in the expansion of dg/dt and dh/dt, and we see that they are identical. 

We can use Equations 196 and 197 for high inclinations without any reservation about 
the convergence relative to y or x Q > if l 0 is not in the neighborhood of the critical values I * 39°, 
or I % 141°. However, a few degrees of deviation of I from 39° or 141° is enough to warrant the 
development of the perturbations of the elements into trigonometric series. 

In the close neighborhood of the critical inclination the theory of resonance must be applied. 
The critical inclinations evidently cannot be obtained easily from the developments of Delaunay. 

It is of interest to mention that in the systems of Jupiter and Saturn the inclinations of the 
orbit planes of the satellites with respect to the orbit planes of the planets all lie between the 
limits given above. The inclination of the orbit plane of Jupiter VUI lies near the resonance case. 

CONCLUSIONS 

The Krylov -Bo go lubov method for eliminating the short-period terms can be applied to solve 
the stellar three-body problem and, in particular, to solve the lunar problem. The perturbations 
in the elements e, c, and l can be obtained in the form of Fourier series with the arguments linear 
in the independent variable. 

The method of Krylov -Bo go lubov can serve to amend the theory of Delaunay because it 
leads to series in a more compact form and their convergence can be investigated more easily. 

In addition, they are applicable to any orbital inclination up to the critical value. 

Our next goal consists in the further application of the Krylov- Bogolubov method to establish 
the solution including the terms of higher order in \ and a. The work will be extended in two di- 
rections: to find the compact form of the analytical development and to arrange the computational 
scheme in order to use the method of successive approximation. The development will be based 
on the author's scheme (Reference 15) of computing the higher order effects in the Krylov- 
Bogolubov method. 


(Manuscript received February 10, 1966) 
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APPENDIX A 


Symbols 

'i - mean anomaly of B, the close companion in ternary system, 
e - osculating eccentricity of orbit of b. 
a - osculating semimajor axis of B. 

7 t - osculating longitude of the perigee of B. 

fi - osculating longitude of the ascending node of B. 

r - position vector of b relative to the central body A, r = | r| . 

e - osculating Laplacian vector of B. 

c - osculating areal velocity of B, c l9 c 2 , c 3 - the notations for the components of c. At the end 
of the exposition the notations and c 2 designate the mean motions of the perigee and of the 
node. 

Cj - mean motion of longitude of the pericenter of B. 
c 2 - mean motion of the longitude of the ascending node of b. 

I - mutual orbital inclination of the close and of the distant companion c. 
r ' - position vector of c relative to the center of masses of A and B. 
f ' - true anomaly of c . 

c' - eccentricity of the orbit of the distant companion C. 

R - disturbing function. 
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